In [ ]:
using LinearAlgebra
using Symbolics
using Plots
using LaTeXStrings
using Distributions
using MAT

Modelo SIR¶

Modelo SIR¶

In [ ]:
N=1
x,y,z,v_x,v_y,v_z,r0,rs,vs,m,t=symN(N)
γ=1/2
β=1/2
conds=[0 0 0 997 3 0]
NN=sum(conds)
F=[-β*v_x[1]*v_y[1]/NN β*v_x[1]*v_y[1]/NN-γ*v_y[1] γ*v_y[1]]

tf=100
dt=0.01

xs,ys,zs,vxs,vys,vzs,ts=RK5(conds,dt,tf,F,false,vs);

plot(ts,vxs',label="Susceptoibles " *L"S_0="*string(vxs[1]),color="Blue")
plot!(ts,vys',label="Infectados "   *L"I_0="*string(vys[1]),color="Red")
plot!(ts,vzs',label="recuperados "   *L"R_0="*string(vzs[1]),color="green")
title!("Solucion Modelo SIR "*L" \beta= "*string(β)*" "*L"\gamma="*string(γ) )
Out[ ]:
No description has been provided for this image
In [ ]:
γ=1/20
β=1/2
conds=[0 0 0 990 10 0]
NN=sum(conds)
F=[-β*v_x[1]*v_y[1]/NN β*v_x[1]*v_y[1]/NN-γ*v_y[1] γ*v_y[1]]

tf=100
dt=0.01

xs,ys,zs,vxs,vys,vzs,ts=RK5(conds,dt,tf,F,false,vs);

plot(ts,vxs',label="Susceptoibles " *L"S_0="*string(vxs[1]),color="Blue")
plot!(ts,vys',label="Infectados "   *L"I_0="*string(vys[1]),color="Red")
plot!(ts,vzs',label="recuperados "   *L"R_0="*string(vzs[1]),color="green")
title!("Solucion Modelo SIR "*L" \beta= "*string(β)*" "*L"\gamma="*string(γ) )
Out[ ]:
No description has been provided for this image
In [ ]:
β=0.345
γ=0.0743#0.157
conds=[0 0 0 5177 23 0]
NN=sum(conds)
F=[-β*v_x[1]*v_y[1]/NN β*v_x[1]*v_y[1]/NN-γ*v_y[1] γ*v_y[1]]

tf=100
dt=0.01

xs,ys,zs,vxs,vys,vzs,ts=RK5(conds,dt,tf,F,false,vs);
plot(ts,vxs',label="Susceptoibles " *L"S_0="*string(vxs[1]),color="Blue",size=(1000,800))
plot!(ts,vys',label="Infectados "   *L"I_0="*string(vys[1]),color="Red")
plot!(ts,vzs',label="recuperados "   *L"R_0="*string(vzs[1]),color="green")
title!("Solucion Modelo SIR "*L" \beta= "*string(β)*" "*L"\gamma="*string(γ) )
Out[ ]:
No description has been provided for this image
In [ ]:
Sr,Ir,Rr=ruido(vxs',vys',vzs',0.0);
# Create a sample matrix
#data = Dict("S" => Sr, "I" => Ir, "R" => Rr, "t" => ts,"N"=>NN)
data = Dict("S" => Sr[1,1:100:end], "I" => Ir[1:100:end], "R" => Rr[1:100:end], "t" => ts[1:100:end],"N"=>NN)

# Specify the file name (change this to your desired file name)
filename = "SIR_limpio.mat"

# Save the data to the .mat file
MAT.matwrite(filename, data)
In [ ]:
Sr,Ir,Rr=ruido(vxs',vys',vzs',0.05);
# Create a sample matrix
#data = Dict("S" => Sr, "I" => Ir, "R" => Rr, "t" => ts,"N"=>NN)
data = Dict("S" => Sr[1:100:end], "I" => Ir[1:100:end], "R" => Rr[1:100:end], "t" => ts[1:100:end],"N"=>NN)

# Specify the file name (change this to your desired file name)
filename = "SIR_ejemp2.mat"

# Save the data to the .mat file
MAT.matwrite(filename, data)
In [ ]:
Is[1]

22.86423174304127
In [ ]:
data = matread("C:\\Users\\Arif\\OneDrive - Instituto Tecnologico y de Estudios Superiores de Monterrey\\Estocasticos\\Reto\\Fase 3\\SIRs1.mat")
Ss=data["Ss"];
Is=data["Is"];
Rs=data["Rs"];
tss=data["ts"];
plot!(tss,Ss',label="Susceptoibles Suavizados " *L"S_0="*string(Ss[1]),color="Blue",line = (:dash, 2))
plot!(tss,Is',label="Infectados Suavizados "   *L"I_0="*string(Is[1]),color="Red",line = (:dash, 2))
plot!(tss,Rs',label="recuperados Suavizados"   *L"R_0="*string(Rs[1]),color="green",line = (:dash, 2))
Out[ ]:
No description has been provided for this image

Pregunta 4 Fase 3¶

In [ ]:
β=0.5785#0.0079
γ=0.0597#0.0008#0.157
conds=[0 0 0 1223 7 0]
NN=sum(conds)
F=[-β*v_x[1]*v_y[1]/NN β*v_x[1]*v_y[1]/NN-γ*v_y[1] γ*v_y[1]]

tf=50
dt=0.01

xs,ys,zs,vxs,vys,vzs,ts=RK5(conds,dt,tf,F,false,vs);
plot(ts,vxs',label="Susceptoibles " *L"S_0="*string(vxs[1]),color="Blue",size=(1000,800))
plot!(ts,vys',label="Infectados "   *L"I_0="*string(vys[1]),color="Red")
plot!(ts,vzs',label="recuperados "   *L"R_0="*string(vzs[1]),color="green")
title!("Solucion Modelo SIR "*L" \beta= "*string(β)*" "*L"\gamma="*string(γ) )
Out[ ]:
No description has been provided for this image
In [ ]:
Is[1]
Out[ ]:
9.040857028704194
In [ ]:
filename = "C:\\Users\\Arif\\OneDrive - Instituto Tecnologico y de Estudios Superiores de Monterrey\\Estocasticos\\Reto\\Fase 3\\matriz_combinada.mat"
mat_data = matread(filename)
Ss = mat_data["matriz_combinada"][1, :]
Is = mat_data["matriz_combinada"][2, :]
Rs = mat_data["matriz_combinada"][3, :]
tss = mat_data["matriz_combinada"][4, :]
plot!(tss,Ss,label="Susceptoibles Suavizados " *L"S_0="*string(Ss[1]),color="Blue",line = (:dash, 2))
plot!(tss,Is,label="Infectados Suavizados "   *L"I_0="*string(Is[1]),color="Red",line = (:dash, 2))
plot!(tss,Rs,label="recuperados Suavizados"   *L"R_0="*string(Rs[1]),color="green",line = (:dash, 2))
Out[ ]:
No description has been provided for this image
In [ ]:
maximum(tss)

100.0
In [ ]:
plot(tss[1:end-1],diff(Ss,dims=2)',label="Suavizado")
plot!(tss[1:end],-(0.345.*(Is.*Ss.*diff(tss,dims=1)[1])./NN)',label="optimizado")
In [ ]:
(Is.*Ss.*diff(tss,dims=1)[1])./N

1×100 Matrix{Float64}:
 1.18372e5  1.23605e5  1.3212e5  …  16261.0  15651.8  15071.7  14531.3
In [ ]:
γ=0.1
β=0.2
conds=[0 0 0 998 1 0]
NN=sum(conds)
F=[-β*v_x[1]*v_y[1]/NN β*v_x[1]*v_y[1]/NN-γ*v_y[1] γ*v_y[1]]

tf=100
dt=0.01

xs,ys,zs,vxs,vys,vzs,ts=RK5(conds,dt,tf,F,false,vs);

plot(ts,vxs',label="Susceptoibles " *L"S_0="*string(vxs[1]),color="Blue")
plot!(ts,vys',label="Infectados "   *L"I_0="*string(vys[1]),color="Red")
plot!(ts,vzs',label="recuperados "   *L"R_0="*string(vzs[1]),color="green")
title!("Solucion Modelo SIR "*L" \beta= "*string(β)*" "*L"\gamma="*string(γ) )
Out[ ]:
No description has been provided for this image
In [ ]:
F'
Out[ ]:
$$ \begin{equation} \left[ \begin{array}{c} - 0.0002002 v_{x_1} v_{y_1} \\ - 0.1 v_{y_1} + 0.0002002 v_{x_1} v_{y_1} \\ 0.1 v_{y_1} \\ \end{array} \right] \end{equation} $$

Plot N -parts¶

In [ ]:
function plotN(xs,ys,zs,N)
gr(size=(800,800))
    if length(xs[1,:])>=1000
        xs=xs[:,1:100:end]
        ys=ys[:,1:100:end]
        zs=zs[:,1:100:end]
    end
        #plot_layout = @layout [a,b]
p = plot()#plot(layout=plot_layout)
if length(size(xs))>1
    N=size(xs)[1]
else
    N=1
end

if N==1
    scatter!([xs[N,end]],[ys[N,end]],[zs[N,end]],label="Posición Final particula "*string(N),color="red",markersize=4,aspect_ratio= :equal,xlabel=L"X \mu m",ylabel=L"Y \mu m")
    plot!(xs[:],ys[:],zs[:],label=false,color="red",lw=2,alpha=0.7) #
else
    for i in 1:N
        #if i>1
        #plot(xs[i,:],ys[i,:],label=false,color=i,lw=2,xlim = [minimum(xs), maximum(xs)],ylim = [minimum(ys), maximum(ys)])
        #end
        plot!(xs[i,:],ys[i,:],zs[i,:],label=false,color=i,lw=2,xlim = [minimum(xs), maximum(xs)],ylim = [minimum(ys), maximum(ys)],xlabel=L"X \mu m",ylabel=L"Y \mu m",aspect_ratio=:equal,alpha=0.7) 
        scatter!([xs[i,end]],[ys[i,end]],[zs[i,end]],label="Posición Final particula "*string(i),color=i,markersize=4)
    end
end
display(p)
end
Out[ ]:
plotN (generic function with 1 method)

Runge Kutta 5to Orden¶

In [ ]:
function RK5(c0,dt,tf,F1,R,vs)
    F1=hcat(vs,F1,zeros(N))
    F=build_function(F1,[x...]...,[y...]...,[z...]...,[v_x...]...,[v_y...]...,[v_z...]...,t,expression=Val{false})
    ##vectores
    ts=collect(0:dt:tf)
    n=length(ts)
    ####################################Posiciones
    xs=zeros(N,n)
    ys=zeros(N,n)
    zs=zeros(N,n)
    
    vxs=zeros(N,n)
    vys=zeros(N,n)
    vzs=zeros(N,n)
    
    #####################################Velocidades
    
#############Condiciones iniciales
    xs[:,1]=c0[:,1]
    ys[:,1]=c0[:,2]
    zs[:,1]=c0[:,3]
    vxs[:,1]=c0[:,4]
    vys[:,1]=c0[:,5]
    vzs[:,1]=c0[:,6]



for i in 1:n-1
        #Runge Kutta
        z0=hcat(xs[:,i],ys[:,i],zs[:,i],vxs[:,i],vys[:,i],vzs[:,i],ones(N).*ts[i])
        
        k1=F[1](z0...)
        k2=F[1](z0.+(dt*k1)/2...)
        k3=F[1](z0.+dt*(3*k1+k2)/16...)
        k4=F[1](z0.+dt*k3/2...)
        k5=F[1](z0.+dt*(-3*k2+6*k3+9*k4)/16...)
        k6=F[1](z0.+dt*(k1+4*k2+6*k3-12*k4+8*k5)/7...)
        
        k1=7*(k1+k6)
        k3=16*k3
        k2=16*k5
        k4=12*k4
        
        xs[:,i+1]=xs[:,i]+(dt/90)*(k1[:,1]+2*k2[:,1]+2*k3[:,1]+k4[:,1])
        ys[:,i+1]=ys[:,i]+(dt/90)*(k1[:,2]+2*k2[:,2]+2*k3[:,2]+k4[:,2])
        zs[:,i+1]=zs[:,i]+(dt/90)*(k1[:,3]+2*k2[:,3]+2*k3[:,3]+k4[:,3])
        
        vxs[:,i+1]=vxs[:,i]+(dt/90)*(k1[:,4]+2*k2[:,4]+2*k3[:,4]+k4[:,4])
        vys[:,i+1]=vys[:,i]+(dt/90)*(k1[:,5]+2*k2[:,5]+2*k3[:,5]+k4[:,5])
        vzs[:,i+1]=vzs[:,i]+(dt/90)*(k1[:,6]+2*k2[:,6]+2*k3[:,6]+k4[:,6])
        
         
        #axs3[i]=h1(xs1[i],xs2[i],xs3[i],ys1[i],ys2[i],ys3[i],zs1[i],zs2[i],zs3[i],vxs1[i],vxs2[i],vxs3[i],vys1[i],vys2[i],vys3[i],vzs1[i],vzs2[i],vzs3[i],ts[i])    
        #ays3[i]=h2(xs1[i],xs2[i],xs3[i],ys1[i],ys2[i],ys3[i],zs1[i],zs2[i],zs3[i],vxs1[i],vxs2[i],vxs3[i],vys1[i],vys2[i],vys3[i],vzs1[i],vzs2[i],vzs3[i],ts[i])
        #azs3[i]=h3(xs1[i],xs2[i],xs3[i],ys1[i],ys2[i],ys3[i],zs1[i],zs2[i],zs3[i],vxs1[i],vxs2[i],vxs3[i],vys1[i],vys2[i],vys3[i],vzs1[i],vzs2[i],vzs3[i],ts[i]);
    end
    
    if R==true
        nn=Int(round(n/1000))
        gr(size=(800,500))
        plot(xs[1,1:nn:end],ys[1,1:nn:end],zs[1,1:nn:end],label="masa 1",color="yellow",lw=5)
        #plot!(xs2[1:nn:end],ys2[1:nn:end],zs2[1:nn:end],label="masa 2",color="blue",lw=2)
        #plot!(xs3[1:nn:end],ys3[1:nn:end],zs3[1:nn:end],label="masa 3",color="red",lw=2)       
    end       

    return xs,ys,zs,vxs,vys,vzs,ts
    ########Graficas####################################
end
Out[ ]:
RK5 (generic function with 1 method)

Animacion 3D¶

In [ ]:
#xs1=
n=length(ts)
mm=5
#xx=vcat(xs1,xs2,xs3)

#yy=vcat(ys1,ys2,ys3)
#zz=vcat(zs1,zs2,zs3)
ext=1/3
l=1 #rotacion de la grafica
anim3d=@animate for i in 1:10:n;
    gr(size=(1000,800))
    a=1
    j=1
    if i>1
        j=Int(ceil(i*0.03))
        a=range(0,1,length=(i-j+1))
        
    end 
    
    Plots.plot(xs[1,j:i],ys[1,j:i],zs[1,j:i],color="yellow",lw=2
    #Plots.plot(xs1[1:i],ys1[1:i],zs1[1:i],color="yellow",lw=2,xlim = [minimum(xx)-ext, maximum(xx)+ext],
    ,label="",xlabel="x",ylabel="y",zlabel="z",camera = (i/4*l, 20)
    ,bgcolor="black",grid=false,alpha=a);
    Plots.scatter!([xs[1,i]],[ys[1,i]],[zs[1,i]], color = "yellow", markersize = mm,label="Masa1");
    
    #Plots.plot!(xs[2,j:i],ys[2,j:i],zs[2,j:i],color="blue",lw=2,label="",alpha=a);
    #Plots.scatter!([xs[2,i]],[ys[2,i]],[zs[2,i]], color = "blue", markersize = mm ,label="Masa2");
    
    #Plots.plot!(xs[3,j:i],ys[3,j:i],zs[3,j:i],color="red",lw=2,label="",alpha=a);
    #Plots.scatter!([xs[3,i]],[ys[3,i]],[zs[3,i]], color = "red", markersize = mm,label="Masa3",markershape=:star4);
    
end;
In [ ]:
gif(anim3d,"jimy .gif",fps=24)

[ Info: Saved animation to C:\Users\Arif\OneDrive - Instituto Tecnologico y de Estudios Superiores de Monterrey\Escuela\Carrera\4to Semestre\mecánica clásica\Reto\jimy .gif
Out[ ]:
No description has been provided for this image

Animacion 2D¶

In [ ]:
n=length(ts)
mm=5
xx=vcat(xs1,xs2,xs3)
yy=vcat(ys1,ys2,ys3)
sc=1/100
anim2d=@animate for i in 1:10:n;
    gr(size=(1000,1000))
    Legend=(bgcolor=(:dark))
    a=1
    if i>1
        a=range(0,1,length=i)
    end 
    
    Plots.plot(xs1[1:i],ys1[1:i],color="yellow",lw=2,xlim = [minimum(xx)-1, maximum(xx)+1],
    ylim = [minimum(yy)-1, maximum(yy)+1],label="",xlabel="x",ylabel="y",bgcolor="black",grid=false,axis=false,alpha=a);
    Plots.scatter!([xs1[i]],[ys1[i]], color = "yellow", markersize = mm,label="Masa1");
    Plots.quiver!([xs1[i]],[ys1[i]],quiver=([axs1[i].*sc],[ays1[i].*sc]),color="yellow",lw=3,length=1/80)
    
    Plots.plot!(xs2[1:i],ys2[1:i],color="blue",lw=2,label="",alpha=a);
    #Plots.scatter!([xs2[i]],[ys2[i]], color = "blue", markersize = mm ,label="Masa2");
    #Plots.quiver!([xs2[i]],[ys2[i]],quiver=([axs2[i].*sc],[ays2[i].*sc]),color="blue",lw=3,length=1/80)
    
    Plots.plot!(xs3[1:i],ys3[1:i],color="red",lw=2,label="",alpha=a);
    #Plots.scatter!([xs3[i]],[ys3[i]], color = "red", markersize = mm,label="Masa3");
    #Plots.quiver!([xs3[i]],[ys3[i]],quiver=([axs3[i].*sc],[ays3[i].*sc]),color="red",lw=3,length=1/80,normalize=true)
    #
end;
In [ ]:
gif(anim2d,"Ciclotron2d.gif",fps=24)

[ Info: Saved animation to C:\Users\Arif\OneDrive - Instituto Tecnologico y de Estudios Superiores de Monterrey\Escuela\Carrera\4to Semestre\mecánica clásica\Reto\Ciclotron2d.gif
Out[ ]:
No description has been provided for this image

Animación¶

In [ ]:
function animac(xs,ys,zs,N,n)
    anim2d=@animate for i in 1:50:n;
        gr(size=(1000,1000))
        a=1
        j=1
        if i>1
            j=Int(ceil(i*0.3))
            a=range(0,1,length=(i-j+1))
            
        end 
        for d in 1:N
            if d==1
                Plots.plot(xs[d,j:i],ys[d,j:i],zs[d,j:i],color=d,lw=2,label="",xlabel=L"X \mu m",ylabel=L"Y \mu m",bgcolor="black",alpha=a,xlim = [minimum(xs), maximum(xs)],ylim = [minimum(ys), maximum(ys)],aspect_ratio=:equal);
                Plots.scatter!([xs[d,i]],[ys[d,i]],[zs[d,i]], color = d,markersize = 6,label="Particula"*string(d),markershape=:star4);
            else
                Plots.plot!(xs[d,j:i],ys[d,j:i],zs[d,j:i],color=d,lw=2,label="",xlabel=L"X \mu m",ylabel=L"Y \mu m",bgcolor="black",alpha=a,xlim = [minimum(xs), maximum(xs)],ylim = [minimum(ys), maximum(ys)],aspect_ratio=:equal);
                Plots.scatter!([xs[d,i]],[ys[d,i]],[zs[d,i]], color = d,markersize = 6,label="Particula"*string(d));
            end
        end
        
        
    end;
    gif(anim2d,"n-cuerpos.gif",fps=24)
end
Out[ ]:
animac (generic function with 1 method)
In [ ]:
function symN(N)
    @variables X,Y,Z,M,x[1:N],y[1:N],z[1:N],v_x[1:N],v_y[1:N],v_z[1:N],t,m[1:N],a_x[1:N],a_y[1:N],a_z[1:N],r
r0=[X,Y,Z]
rs=[[x...]';[y...]';[z...]']'
vs=[[v_x...]';[v_y...]';[v_z...]']'
    return x,y,z,v_x,v_y,v_z,r0,rs,vs,m,t
end
Out[ ]:
symN (generic function with 1 method)
In [ ]:
function ruido(S,I,R,p)
    N=length(S)
    ch=[1,0,-1]
    for i in 1:N
        r1=rand(ch,4).*p
        S[i]=S[i]+r1[1]*(S[i])
        #E[i]=E[i]+r1[2]*(E[i])
        I[i]=I[i]+r1[3]*(I[i])
        R[i]=R[i]+r1[4]*(R[i])
    end
    return S,I,R
end

ruido (generic function with 1 method)